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Abstract. We explore the ability of directional nuclear-recoil detectors to constrain the 
local velocity distribution of weakly interacting massive particle (WIMP) dark matter by 
performing Bayesian parameter estimation on simulated recoil-event data sets. We discuss 
in detail how directional information, when combined with measurements of the recoil- 
energy spectrum, helps break degeneracies in the velocity-distribution parameters. We 
also consider the possibility that velocity structures such as cold tidal streams or a dark 
disk may also be present in addition to the Galactic halo. Assuming a CF4 detector 
with a 30-kg-yr exposure, a 50-GeV WIMP mass, and a WIMP-nucleon spin-dependent 
cross-section of 10~^ pb, we show that the properties of a cold tidal stream may be well 
constrained. However, measurement of the parameters of a dark-disk component with a 
low lag speed of ^^50 km/s may be challenging unless energy thresholds are improved. 
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1 Introduction 

Solid-state and liquid WIMP-dark-matter detectors designed to measure the energy of 
nuclear recoils from WIMP collisions are entering maturity on both theoretical and experi- 
mental fronts. Indeed, a large number of theoretical studies have investigated the statistical 
power of these experiments to characterize WIMP dark matter [1-8]. Furthermore, a va- 
riety of experiments [9-19] are currently running, with a few tantalizing signals already 
observed [20-22]. In contrast, gas detectors with sensitivity to the nuclear-recoil direction 
via the measurement of ionization tracks are still relatively nascent. Nevertheless, some 
theoretical studies on directional dark-matter detection have likewise been conducted [23- 
38], and a small number of directional detectors are currently under development [39-44]. 

A primary advantage of these directional detectors is that they allow the possibility of 
easily distinguishing between terrestial background events (which should be isotropic) and 
WIMP-induced recoil events (which should be non-isotropic, due to our motion through the 
Galactic halo). Perhaps even more intriguing is the possibility that directional detectors 
may allow the details of the local WIMP velocity distribution to be inferred. Theoretical 
expectations and N-body simulations of the Galactic halo both give us reason to believe 
that even if the local spatial distribution of dark matter might be expected to be rel- 
atively smooth, the velocity distribution may possess interesting structure. Possibilities 
include cold tidal streams passing through the local solar neighborhood, a dark-matter 
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disk aligned with the steUar disk, and warm debris flows [45-60] . Directional-detection ex- 
periments may confirm this picture, perhaps shedding light on not only the local velocity 
distribution but also the process of structure formation on galactic scales. Furthermore, 
a better understanding of the dark-matter velocity distribution will yield improved con- 
straints on the particle properties of the dark matter - in particular, the WIMP mass and 
the WIMP-nucleon cross section. 

In this paper, we explore the statistical ability of directional detectors to constrain the 
local WIMP velocity distribution. The organization of the paper is as follows. In section 2, 
we briefly review the general formalism for obtaining the directional recoil spectrum from 
the dark-matter velocity distribution. We then discuss how a binned likelihood analysis of 
the directional recoil-event data may be used to estimate the parameters of the velocity 
distribution. In section 3, we perform parameter estimation on simulated data sets to 
demonstrate the power of these methods. We consider three specific distributions as exam- 
ples: (1) the standard halo model, (2) a halo model with an additional cold dark-matter 
stream component, and (3) a halo model with an additional dark-matter disk component. 
We discuss implications for future directional dark-matter-detection studies and give our 
conclusions in section 4. 



2 Formalism 

2.1 The Directional Recoil Spectrum 

The directional dependence of the WIMP signal was first discussed in ref. [61]. The for- 
malism presented here for the calculation of the nuclear-recoil event rate (as a function of 
recoil energy and direction) induced by WIMPs from a given WIMP velocity distribution 
was worked out in ref. [62]. 

Let the distribution of WIMP velocities Vg in the Galactic rest frame be given by 
/g(vg). If we neglect the gravitational influence of solar-system bodies, the distribution of 
WIMP velocities v in the lab frame moving with velocity viab with respect to the Galactic 
frame is then 

/(v) = /g(v + vi,b) , (2.1) 

since the various velocities are related by v = Vg — viab- 

It can be shown that the directional recoil rate is given by 

dR _/>o--5(^);(,^,g). (2.2) 



Here, R is the number of events per exposure (detector mass multiplied by time), po 
is the local WIMP density, fi]\f = m^m^/{m^ + m^) is the reduced mass of the WIMP- 
nucleus system, is the lab-frame nuclear-recoil momentum, E = g^/2mN is the lab-frame 
nuclear-recoil energy, = q/2ii is the minimum lab- frame WIMP speed required to yield a 
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recoil energy E, the WIMP-nucleus elastic cross section is da/dq^ = a]\fS{q)/4i2'j^v'^ (where 
S{q) is the nuclear form factor), and 

/K,q) = |5{^r■n-v^)f{^r)d'v (2.3) 

is the Radon transform of the lab-frame WIMP velocity distribution. In practice, it is 
often easier to take the Radon transform of the Galactic-frame velocity distribution, and 
then use the relation 

f{Vci, q) = fg{Vq + Vlab " Q, Q) (2.4) 

in eq. (2.2) to find the directional recoil spectrum. In this paper, we shall assume that 
the velocity of the lab frame is given by that of the local standard of rest (LSR); i.e., we 
take Vlab = (^LSRi ^LSRi ^lsr) = (220 km/s,90°,0°) in Galactic coordinates [63]. We shall 
ignore both the motion of the Earth around the Sun and the rotation of the Earth. 



2.2 The Binned Likelihood Function 

Our ultimate goal will be to investigate the degree to which a binned likelihood analysis 
of observed recoil events might recover the parameters of the dark-matter velocity distri- 
bution. Given a number of observed events, we may construct a sky map of the data by 
binning both signal and background events into A'pix pixels (we shall use the pixels of equal 
angular area given by HEALPix [64] in this paper), as well as binning events in each pixel 
into A^bins energy bins (we shall also assume the energy bins are of equal width). That is, 
the sky map specifies the number of observed signal and background events Mjj in the i-th 
pixel of solid angle dQi and the j-th energy bin of width AEj for all A'^pix pixels. We use 
pixels and bins as a proxy for finite angular and energy resolution in the detectors. 

We would then like to construct a likelihood function that may be used to compare 
the observed sky map to the predicted sky map. The predicted sky map can be specified 
by the total number of events A'tot) which are observed over the energy-sensitivity range 
[Erain, ErasLx], as Well as the normalized distribution P{q,E) in angle and energy of these 
events. These quantities are simply related to the direction recoil spectrum via 

dn 

XNtotPiq,E)=£—-—, EG 
dEdilq 

Here, A is related to the background-rejection power of the experiment, defined so that 
the number of signal events is A'^gig = AA'tot and the number of background events is 
A^'bg = (1 — A)A'^tot- The effective exposure £ gives the fraction of the total exposure 
£tot arising from the total mass of target nuclei; we have implicitly assumed the detector 
acceptance is not energy dependent. 
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Note that the predicted angular distribution of signal events is then given by a nor- 
malized integral over the energy-sensitivity range 



P(q) 



(2.6) 



where Emm = '^v'^,mmf^N / '^^ and E^max = 2vl -m.^^fi%/mt^. The energy distribution P{E) 
of signal events is similarly given by a normalized integral over all angles. 

For a given velocity distribution, we see that the predicted mean number of signal 
and isotropic background events in the i-th pixel (centered at the direction qj) and j-th. 
energy bin is then given by 



where Pb{E) is the energy distribution of the isotropic background events, normalized to 
unity over the recoil-energy sensitivity range. The approximation in the second line (which 
approximates the angular integral over each pixel with the value of the integrand at the 
center of each pixel multiplied by the pixel size) holds in the limit that A'pix is large. In each 
energy-binned pixel, the number of events is Poisson distributed, so a suitable likelihood 
function is given by the product of the distributions in each energy-binned pixel 



where p{Mij\Mij) is the Poisson distribution function for the random variable Mij (the 
observed number of events) with mean Mij (the predicted number of events, which depends 
on the velocity distribution /g and the experimental conditions). 

Let us now assume that the Galactic-frame velocity distribution /g depends on the 
model parameters 0^, the measurement of which is of interest. Then the Radon transform 
/(vq,q) of the lab-frame velocity distribution also depends on 0^. From eq. (2.4), we see 
it further depends on viab, which may be treated as three additional model parameters 
{'t'lab) ^lab) ^lab} (in Galactic coordinates). Thus, it is clear that P{^,E) depends on 0^ and 
viab- Furthermore, examining the form of P(q, i?) in eq. (2.6), we see that it additionally 
depends on /xat through the argument of the form factor and the limits of the integrals. 
Therefore, it depends on the known target-nucleus mass ttt-n and the unknown WIMP mass 




(2.7) 



(2.8) 



i=i j=i 
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rriy^, the latter of which may be treated as yet another model parameter we would like to 
measure. However, the direct dependence of i-*(q, E) on weakens at high values of m^, 
since the dependence of iin on is also weak at such values. Finally, the background- 
rejection power of the detector may not be well known, so we may also treat A as a model 
parameter. 

The likelihood function C can therefore be written in terms of the functions and 
parameters specifying both the velocity distribution and the experimental conditions as 
C[m^,X,via.h,ha.hAi>.h,dk,S{q); mN,Eram,Eraa.^,Npi^,Nuns,Ntot,PB{E)]. By assuming the forms 
of S{q) and Pb{E), a likelihood analysis of observed data can then be done to estimate 
the parameters that are unknown. In this work, we shall take S{q) = 1 and assume that 
Pb{E) = l/(£'max — Erain) IS a flat spectrum for simplicity. 

Note that in this treatment, we have folded all the dependence of the amplitude of 
the signal A'^sig = AA'tot into A. However, from eqs. (2.2) and (2.5), it is clear that A'^gig will 
itself depend not only on the exposure £tot, but also on the additional, totally degenerate 
parameters po and a^. If one is more interested in the particle properties of the WIMP, 
the usual procedure is to fix pQ « 0.3 GeV/cm^ (a typical estimate of the local dark- 
matter density, which may be measured by various means [65-69]) and to treat a at as the 
parameter of interest. In fact, more often the focus is placed on the WIMP-nucleon cross 
section ap^n', for a spin-dependent interaction, this is related to the WIMP-nucleus cross 
section via 

CTN _ 4 /i^ J + 1 {ap{Sp) + an{Sp)f 

fj ~ 3 /;2 / o2 ' y'^-^> 

where fip^n is the WIMP-nucleon reduced mass, J is the nuclear spin, ap^n are the effective 
nucleon coupling strengths, and {Sp^n) are the expectation values of the spin content of the 
nucleon group [70, 71]. 

However, it can be argued that the typical value of po ^ 0.3 GeV/cm"^ usually assumed 
may not even be relevant for direction-detection experiments; it is simply a large-scale av- 
erage of the dark-matter density at the Galactic radius of the Sun ro ~ 8.5 kpc, and does 
not account for the possible existence of substructure in the immediate neighborhood of the 
Earth [50, 72]. The assumption of this typical value then strongly colors any conclusions 
drawn about the estimated value of (Tp^„. Thus, perhaps a more assumption-independent 
approach would be to ask: if we have observed a given number of events A'totj how well 
can the parameters of an assumed velocity distribution be estimated using only the distri- 
bution P(q, E) of these events? In this way, we will sidestep the issues introduced by the 
degeneracy of po and a^, which are now subsumed into the single parameter A. 

3 Likelihood Analyses of Simulated Data 

We shall now explicitly demonstrate the feasibility of estimating the parameters of the 
velocity distribution by performing likelihood analyses on simulated data sets. We shall 
consider three parameterized velocity distributions: 1) the standard halo model, 2) a halo 
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model with an additional cold stream component, and 3) a halo model with a disk com- 
ponent. The method of analysis is as follows. Assigning fiducial values for the velocity- 
distribution and experimental parameters, we first randomly generate a number of recoil 
events using the procedure described in appendix A. We then bin the simulated events in 
angle (using HEALPix) and in energy to create a simulated sky map. We then use this to 
calculate the likelihood function, employing MultiNest [73, 74] to sample the likelihood 
function within the model parameter space, assuming flat priors. The GETDIST routine 
from the COSMOMC package [75] is then used to calculate the ID and 2D marginalized 
posterior probability distributions. We also calculate the minimum credible intevals (MCIs) 
(as defined in [76]) of the posterior probability distributions, and examine how well the 
fiducial parameters have been recovered. 

The fiducial values chosen for the velocity-distribution parameters will be discussed 
below for each of the three cases. However, let us first motivate the choice of the values for 
the experimental parameters. We shall perform the analyses assuming that the simulated 
data were collected by a CF4 MIMAC-like experiment [39, 77-79]. That is, we assume the 
target nucleus is ^'^F, and that the relevant WIMP-nucleus interaction is spin-dependent 
and can be modeled using J = 1/2, and assuming a pure proton coupling with Op = 1, 
ttn = 0, and (Sp) = 0.5. Furthermore, we take the energy-sensitivity range to be 5- 
50 keV. This range corresponds to that quoted by the MIMAC collaboration; the 5-keV 
threshold arises from the ionization threshold (taking into account quenching), while the 
upper bound is chosen to limit contamination from background events that dominate the 
signal at higher energies. We shall take a fiducial value of = 50 GeV, so we see that 
the 5 keV recoil-energy threshold corresponds to a sensitivity to WIMP velocities down to 
~150 km/s. 

The experimental parameters also include the angular and energy resolution. In all 
of the analyses below in which we use directional information, we shall take the number 
of pixels to be the same, setting A'^pix = 768 (i.e., HEALPix order 8). This roughly 
corresponds to the ~ 10° angular resolution expected to be attained by a MIMAC-like 
experiment. On the other hand, we shall investigate the effect of varying the energy 
resolution by considering different numbers of energy bins A^bins for each analysis, which 
will be given in detail below. For now, we note that the bin widths we shall assume are 
relatively large and conservative, considering that the energy resolution of the micromegas 
detectors used in the MIMAC experiment is expected to be ~15%. In any case, the exact 
values assumed for the angular and energy resolutions do not have a large effect on the 
quality of the parameter estimation, and it can be shown that unbinned likelihood analyses 
yield similar results. 

Finally, there remains the question of the number of signal and background events we 
should examine for each analysis. We will assume various values of Ns[g and Ni,g for each 
case, as will be discussed and motivated below. 

Values for the experimental parameters used for each analysis are summarized in 
table 1. Fiducial values for the velocity-distribution parameters and the flat prior ranges 
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model 


iVpix 


^bins 


l^Kev j 


-E'max 

l^Kev j 


A^sig 


iVbg 


iVtot 


li^s yi ) 


I [Kg-yrj ) 


(m^ fixed) 


7fi8 
(or 1) 


1 n 

(or 1) 


5 


50 


100 





100 


4.4 





lLCil\J \JL11 y 

(m^ flat prior) 


7fi8 
(or 1) 


10 

(or 1) 


5 


50 


100 





100 


4.4 





halo-only 

(6 parameters) 


768 
(or 1) 


10 


5 


50 


100 


100 


200 


4.4 


23 


halo+stream 


768 


20 


5 


50 


650 


300 


950 


28.9 


10.4 


halo+disk 


768 


36 


5 


50 


541 


272 


813 


20 


15 


halo+disk 
(zero threshold) 


768 


40 





50 


900 


300 


1200 


20 


15 



Table 1. Experimental parameters used to simulate data for each analysis. The quoted values of 
£tot and Rhg assume a CF4 detector and the typical values pa 0.3 GeV/cm^ and crp.„ w 10"'^ pb. 



are summarized in table 2; we shall proceed to discuss the choice of these fiducial parameters 
for each velocity distribution in detail. 

3.1 Halo-only Model 

For simplicity, we shall first consider a Galactic dark-matter halo with a velocity distribu- 
tion that may be locally approximated as an isotropic Maxwellian with velocity dispersion 
cth, truncated at the Galactic escape speed Vesc- This truncated-Maxwellian distribution 
is generally referred to as the standard halo model, and is of the form 



/™(Vg;crH,'(;esc) 



1 



iVesc(27ra2)3/2 



exp 



2.^ 



eive 



(3.1) 



where 



erf 



/ Vesc \ 

\V2au) 



2Ve 

vr an 



exp 



24 



The Radon transform of this distribution is given by 



f™{w,w;an,Vc 



/'^(vg 



• W - w)f™{Vg; CTH, Vesc) d^^Vg 



1 



iVesc(27rCT2)l/2 



exp 



w 



exp 



2ct|, 



(3.2) 
(3.3) 

e{vesc - w) (3.4) 



The lab frame moves with respect to this velocity distribution with a velocity viab, so as 
in eq. (2.4) we may find the Radon transform of the velocity distribution in the lab frame 



f^{Vq, q) = fl^iVq + Viab • q, q; O-H, Vesc) , 



(3.5) 
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which yields the directional recoil spectrum via eq. (2.2). 

Note that integration of /^(fq,q) over angles gives the recoil spectrum in the usual 
way; this integral takes the simple analytic form [80] 



dR ^/tt fJH , p 
-TF^ oc — 1= <; erf 

dE 2V2i'iab 



Vq{E) + t;iab 



erf 



V2<TH 



exp(-^), (3.6) 



2cr 



which approaches a falling exponential in the limit that ujab and v^sc ~^ oo. 

For the standard halo model, we see that the set of parameters determining the 
Galactic-frame velocity distribution is simply 9k = {o"h, "fcsc}- Here and afterwards, we shall 
assume Vesc = 550 km/s is known independently (and in practice, the energy-sensitivity 
range is such that the analysis is not sensitive to the exact value of Vgsc)- In total, there are 
then six parameters {m-^. A, uiab, ^iab> ^lab, o"h} that will determine the energy and angular 
distribution of events. In the following analyses, we shall consider a standard halo model 
with fiducial parameter values {m^ = 50 GeV, vi^b = 220 km/s, ^lab = 90°, feiab = 0°, cth = 
viah/V^ = 155 km/s}. Note that the choice of ch = ^lab/v^ yields a velocity distribution 
corresponding to a singular isothermal sphere, with halo profile p{r) w Po{ro/r)'^. 

3.1.1 viah — o"H Analyses 

To first gain a qualitative and intuitive understanding of the additional power that direc- 
tional information provides, let us perform a simple illustrative exercise. Using the fiducial 
parameter values defining the standard halo model above, we generate A'^gig = 100 signal 
and A^bg = background events (i.e., we take A'tot = 100 and A = 1). The spectrum and 
recoil map for these events is shown in figure 1. 

We then further assume that the values of the parameters {m^, A, li^h, hah} are known 
exactly, so that only {fiabj^H} are unknown and remain to be estimated. Sampling the 
likelihood function over the 2D {uiab, ch} parameter space using MultiNest (assuming flat 
priors over the ranges given in table 2), we then perform three separate likelihood analyses 
of these 100 events. For the first analysis, we use only the energy information of the recoil 
events, accomplished by setting A'^pix = 1 and A'^bins = 10. For the second analysis, we use 
only the directional information, setting A'^pix = 768 and A'^bins = 1- Finally, we analyze 
the data using both direction and energy information, setting A'^pix = 768 and A^'bins = 10. 
Doing these three separate analyses allows us to study exactly how the direction and energy 
information translate into information about t^iab and uh- 

The results of the three analyses are presented in the top row of figure 2. We see 
that the energy-only analysis yields contours for the 2D marginalized posterior probability 
distribution in f lab — space that indicate that f lab and an are anti-correlated. This can 
be understood simply by noting from eq. (3.6) that the dependencies of the shape of the 
energy spectrum on vi^h and an are similar. Increasing the value of either viab or cjh results 
in a larger fraction of recoil events at higher energies, with both actions flattening out the 
exponentially- falling recoil spectrum; this behavior is illustrated in flgure 3. Thus, the free 
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parameters vi^h and cjh can only vary in a roughly inverse manner with each other if the 
observed shape of the spectrum is to be maintained. 

On the other hand, the direction-only analysis yields contours that indicate fiab and 
(Th are correlated. Again, this can be easily understood by considering the dependence of 
the directional recoil map on v\ah and cth, illustrated in figure 4. Note that in the limit that 
uiab vanishes, the recoil map becomes isotropic; conversely, it is clear that increasing the 
value of flab makes the map more anisotropic and asymmetric by increasing the enhance- 
ment in one hemisphere due to forward scattering by incoming particles from the "WIMP 
wind." In contrast, the recoil map becomes more isotropic as o"h increases and becomes 
much larger than viab, since the number of incoming WIMPs arriving in the opposite di- 
rection from the WIMP wind is then increased. We see that vi^h and cth must vary in a 
roughly proportional manner with each other if the observed large-scale anisotropy of the 
recoil map is to be maintained. 

Thus, the energy-only and direction-only analyses provide orthogonal sets of informa- 
tion on the velocity and dispersion parameters. It is then easy to see how combining both 
sets of information in the direction+energy analysis yields contours that demonstrate that 
uiab and (Th are relatively uncorrelated. 

We can repeat this exercise using the same recoil-event data set, this time relaxing 
the assumption that the WIMP mass is known independently. We now allow to 
be an additional free parameter to be estimated, assuming a flat prior as shown in table 2 
and sampling over the 3D {m^^-, ?;iab, o"h} parameter space. The results of the energy-only, 
direction-only, and direction+energy analyses are shown in the bottom row of figure 2. 
Interestingly, we see that the quality of the contours in the energy-only analysis is severely 
degraded compared to the case where the mass is known exactly, with a long, flat tail in 
the cth direction appearing in the posterior probability distribution. This can be explained 
in a similar manner as before; allowing for decreased values of is only possible if 
increased values of either vi^h and an compensate to flx the observed fraction of events 
at high energies. However, the slope of the spectrum is slightly more sensitive to changes 
in o"H than to changes in uiab) as can be shown by examining eq. (3.6) (in particular, by 
considering the absolute magnitude of the derivatives of dR/dE with respect to fiab and 
(Th in the relevant regions of parameter space). Thus, small values of are more easily 
compensated for by increasing cjh, so a long tail appears in the (Th direction. 

It is notable that the analyses incorporating directional information are relatively 
insensitive to the lack of prior knowledge of the WIMP mass. This again demonstrates the 
power and robustness of combining directional and energy information to fix parameters 
of the velocity distribution. 

3.1.2 6-parameter Analyses 

We now proceed to perform likelihood analyses over the full 6D standard-halo-model pa- 
rameter space of {m-x' "^labi ^lab; ^labj ('"h}- Using the same fiducial values for the velocity- 
distribution parameters as before, we simulate a total of A'tot = 200 events and take 
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A = 0.5, corresponding to A'^gig = 100 and A'^bg = 100. For a MIMAC-like experiment 
with 10 kg of target mass, assuming a WIMP mass = 50 GeV, a typical local density 
po ~ 0.3 GeV/cm^, and a spin-dependent WIMP-nucleon cross section of ap^n ~ 

10-3 pb 

(consistent with current observational limits from both direct-detection and neutrino exper- 
iments [81-83]), this number of signal events roughly corresponds to a 5- month observation 
period. The number of background events then corresponds to the assumption of a rela- 
tively high background event rate of ~23/kg/yr (in comparison, it may be reasonable to 
expect background event rates as low as 10/kg/yr [36]). We again take A^pix = 768, and 
bin events into A^bins = 10 energy bins. The recoil spectrum and binned recoil sky maps 
for the simulated data are shown in figures 5 and 6, respectively. 

We shall perform both an energy-only (A^pix = 1) analysis and a direction+energy 
analysis. The results of the parameter estimation are shown in the triangle plots of the ID 
and 2D marginalized posterior probability distributions in figures 7 and 8. Note that the 
68% MCIs for each parameter are plotted there, and are also given in table 2 along with 
the ID marginal posterior modes. 

From the triangle plots, it is clear that the analysis incorporating the directional 
information is able to recover the parameters with greater fidelity than the energy-only 
analysis. Not only does this information allow for the direction (Ziabi hab) of the lab frame 
to be recovered quite accurately, it also allows a rough measurement of the background- 
rejection power A from the data itself. That is to say, the directional information indeed 
allows the isotropic background component to be separated from the anisotropic signal 
component. In contrast, the energy-only analysis is unable to recover A correctly. As we 
saw above, flattening of the exponentially-falling spectrum may be caused by variations in 
the other parameters, and cannot be easily separated from the introduction of a truly flat 
background spectral component using only the energy information (at least, not without 
improved statistics from a larger number of events). 

Note also that the estimate of is quite poor for both analyses, with a long tail 
extending to large values of m^^. This partially stems from the fact that we are only us- 
ing information from the distribution of the events, which only has weak dependence on 
as rri), increases, as discussed previously. Since there is some additional dependence 
of the amplitude of the signal on m^, as can be seen from eq. (2.2), making the afore- 
mentioned assumptions about the additional amplitude- fixing parameters pQ and can 
greatly improve these estimates of m^. 

Finally, from examination of the plots for uiab and cth in figures 7 and 8, we see 
that the intuitive results about the ability to constrain these parameters from the simple 
2-parameter and 3-parameter analyses are basically borne out even in the full 6-parameter 
analysis. However, we note that the presence of a flat isotropic background at this level 
does degrade the ability to pin down cjh, even with directional information, as is evident 
from the long tail in the marginalized posterior probability distribution for cth- This tail 
may likewise be reduced by improved statistics or better background rejection. 

To summarize, for the standard halo model, 100 signal events (which might reasonably 
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be observed in a 5-month period) yield a good measurement of the direction of the LSR, and 
allows rough constraints to be placed on the velocity and dispersion of the Galactic halo, 
the presence of a non-negligible background notwithstanding. Of course, the results may 
be improved simply by lengthening the observation period; a 3-year observation period 
resulting in a 30-kg-yr exposure might be expected. In such a period, for the fiducial 
parameter values we have assumed for the standard halo model and WIMP properties, 
one would expect ~650 signal events within the energy-sensitivity range 5-50 keV (and 
~900 total events down to zero threshold). It is then interesting to ask whether or not any 
interesting structure in the local dark-matter velocity distribution, such as cold streams or 
disk components, may be detected with a comparable number of events. We shall proceed 
to investigate this question in the subsections to follow. 

3.2 Halo-l-stream Model 

Motivated by the results of the simulations mentioned in section 1, we next consider a 
model with a local dark-matter stream in addition to the dark-matter halo. We shall 
assume that the velocity distribution in the Galactic frame of a stream component is given 
by 

/g (vg) = /g™(vg - vs; as, ^;osc) , (3.7) 

where the velocity vector of the stream in the Galactic frame is vs (and may be indicated 
by its magnitude vs and its direction (^S)^s) in Galactic coordinates) and the velocity 
dispersion fjs is small for a cold tidal stream. We further assume that the dark-matter 
particles in the stream consist of a fraction of the total particles in both the halo and 
the stream locally. Using the linearity of the Radon transform and eq. (2.4), we find 

/^+^(^^q,q) = (l-^s)/™(^^q + Vlab-q,q;0-H,'i;esc) + ^s/™(?^q + (viab " Vs) • q, q; CJs , Ucsc) , 

(3.8) 

which yields the directional recoil spectrum via eq. (2.2). 

We shall assume the fiducial values of the parameters {m^, fiab, ^lab, ^labj <7h} deter- 
mining the halo component of the distribution are identical to those used in the halo-only 
model above. For simplicity, we shall further assume that the parameter uiab = vi^SR = 
220 km/s is known exactly, so that the halo component is specified by 4 free parame- 
ters. This is done simply to differentiate between the two Maxwellian components of the 
velocity distribution within the likelihood analysis, which would otherwise be arbitrarily 
assigned. By fixing viab) we can identify the "halo" as the Maxwellian component that 
moves with mean velocity fiab with respect to the lab frame, while the "stream" is identi- 
fied as the secondary component. For the 5 stream parameters, we shall take the fiducial 
values {As = 0.1, Is = 65°, 6s = 25°, vs = 510 km/s, as = 10 km/s}. Including A, the 
halo-l-stream model is specified by 10 parameters in total. 

For this halo-|-stream model, we shall perform only a direction-l-energy analysis, as- 
suming we have observed an number of events comparable to that expected in the baseline 
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halo-only scenario with a 30 kg-yr total exposure. For the experimental parameters, we 
adopt {iVpix = 768, iVbins = 20, iVtot = 950, A = 650/950 « 0.684}. This yields iVsig = 650 
and A'^bg = 300, roughly corresponding to an exposure of 28.9-kg-yr and background rate 
of 10.4/kg/yr. These values are close to those in the baseline model, since the 10% stream 
component is merely a small perturbation of the standard halo-only model. 

The recoil spectrum and sky maps for the simulated data are shown in figures 9 and 10, 
respectively. The triangle plot of the posterior probability distributions resulting from the 
direction-l-energy analysis is shown in figure 11, and the posterior modes and MCIs are 
given in table 2. The quality of parameter estimation is quite good; with the exception 
of m^, all the parameters are recovered accurately without bias and with a low degree of 
correlation. 

This seems to suggest that streams may be detectable with a MIMAC-like experiment. 
However, let us note that we have adopted a somewhat unrealistically high stream fraction 
here (solely for illustrative purposes, e.g., to make visible in the recoil maps in figure 10 
the feature arising from the stream). More realistically, simulations suggest that typical 
stream fractions are more likely at the ~1% level, and furthermore, that the probability 
that a tidal stream dominates the local neighborhood of the Sun (averaged on kpc scales) 
is less than 1% [50, 56]. Nevertheless, our analysis here serves as a proof of principle for 
the possibility of detecting structure in the velocity distribution beyond the standard halo 
model and using directional information to constrain parameters. Furthermore, for a more 
realistic stream fraction of Aa, ~ 1%, one might expect only a ~\/lO reduction in statistical 
power when constraining the stream parameters. Finally, this result suggests that one 
would expect comparable - or even improved - statistical power when constraining warm 
debris flows, which may have similarly large mean velocities but are thought to compose 
a more significant fraction (tens of percent, becoming dominant at large velocities) of the 
local dark-matter density. 

3.3 Halo-l-disk model 

Finally, we consider a Galactic model with a dark-matter-disk component in addition to 
the dark-matter halo [49, 51]. We assume that disk rotates such that local particles in 
the disk move with some average velocity vd with respect to the halo/Galactic frame; 
this velocity may be specified by its magnitude vn and its direction (Zd,&d) in Galactic 
coordinates. We shall take vd to be parallel to vlsr, so that these particles lag the LSR 
by via_g = vj) — t'LSR; typical values from simulations are ~ 170 km/s and uiag ~ — 50 
km/s, such that the dark disk rotates more slowly than the stellar disk. We again assume 
that the local velocity distribution of the disk component is also a truncated Maxwellian 
with dispersion cd and that particles in the disk compose a fraction Ad of the local 
particles, so that the disk velocity distribution is identical to that of halo-|-stream model 

for {As,vs,o-s} {Ay),vi),(7b}- 

We again use the fiducial values of the halo component parameters, assuming fiab = 
^^LSR = 220 km/s is known as before. For the disk parameters, we shall take the fiducial 
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values {^D = 0.5, /d = 90°, 6d = 0°, vb = 170 km/s, = 100 km/s}. 

It is clear that detecting nuclear recoils from collisions with low-velocity WIMPs in 
the disk component may be challenging, since the average lab-frame speed = 170 km/s 
for such WIMPs is close to the ~150 km/s velocity threshold (corresponding to the 5 keV 
energy threshold, with our fiducial values for mN and m^). Therefore, to investigate the 
effect of the energy threshold on the parameter estimation, we shall perform two direc- 
tion-l-energy analyses: the first with a zero-energy threshold, and the second with a 5-keV 
threshold as before. 

Again, we would like to evaluate the statistical power of our halo-|-disk analysis, 
assuming we have observed an number of events comparable to that expected in the baseline 
halo-only scenario with a 30-kg-yr exposure. As mentioned previously, this baseline scenario 
would yield ~900 signal events in the 0-50 keV range, along with 300 background events 
(assuming the expected rate of 10/kg/yr). We therefore adopt we adopt {A'^pix = 768, 
A'bins = 40, A'^tot = 1200, A = 0.75} for the first zero-threshold analysis. For the second 
analysis, we then make a cut of all events below the 5-keV threshold, resulting in {Abins = 
36, Atot = 813, A = 541/813 ~ 0.665}. Worth mentioning here is the fact that simulations 
show that the presence of a disk component enhances the local dark-matter density over 
that from the halo alone, so that po ^ po/i^ — ^d)- This correspondingly decreases the 
amount of exposure needed to reach the baseline of 900 signal events, and also allows for 
a larger background rate. 

The recoil spectrum and sky maps are shown in figures 12 and 13, the triangle plots 
for the zero-threshold and 5-keV-threshold analyses are shown in figures 14 and 15, and 
the posterior modes and MCIs are again given in table 2. From the triangle plots, it is 
immediately clear that the finite energy threshold severely restricts the ability to constrain 
the disk-component parameters, even leading to multimodal posterior probabilities for some 
of the parameters. Even in the idealized case of zero threshold the parameter estimation is 
imperfect, as is evidenced by the amount of bias in the estimation of the disk fraction ^d- 
Although these results may be improved with better statistics or reduced background, it 
is clear that measurement of the parameters of a dark-matter disk will be challenging with 
current energy thresholds, if predictions of the disk parameters from N-body simulations 
are indeed valid. 

4 Conclusions 

Using binned likelihood analyses of simulated WIMP-nucleus recoil events, we have demon- 
strated how detectors with directional sensitivity may place constraints on the parameters 
of the local dark-matter velocity distribution. Directional sensitivity allows isotropic back- 
ground events to be distinguished from signal events, but more interestingly, it also helps to 
break degeneracies in the velocity-distribution parameters that cannot be broken with spec- 
tral information alone. As an illustrative example, we examined the case of a Maxwellian 
standard-halo-model velocity distribution. By comparing the statistical power of energy- 
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only, direction-only, and direction+energy analyses, we demonstrated how the degeneracy 
in the standard halo model between the velocity t>iab of the Earth through the halo and 
the halo velocity dispersion uh is broken with directional information. 

If the local dark-matter density is indeed pQ ~ 0.3 GeV/cm^, for a 50 GeV WIMP and 
a WIMP-nucleon spin-dependent cross section of 10"^ pb, a MIMAC-like, 10-kg CF4 di- 
rectional dark-matter detector might observe several hundred WIMP-nucleus recoil events 
in an observation period of 3 years. With 650 signal events, the presence of a cold stream 
composing a fraction of the local dark-matter density may be detected, and its properties 
measured with relatively good accuracy. However, the detection of a dark-disk component 
requires sensitivity to WIMPs moving at relatively low velocities, and may only be feasible 
if energy thresholds are improved. 

We have focused on an approach that seeks to maximize the amount of information 
about the velocity distribution, while minimizing the number of assumptions about the 
particle properties of the WIMP or the local density of dark matter. The only assumptions 
of our approach are: 1) the form of the velocity distribution is known and can be parame- 
terized, 2) the background is flat and isotropic. Future work might consider the relaxation 
of these assumptions. Furthermore, the analyses presented in this paper simply provide 
a proof of principle, showing that the parameter estimation might be feasible in cases of 
interest. However, a more in-depth study of how the statistics of the parameter estimation 
improve as the experimental conditions are varied might be warranted, with the goal of 
finding the best balance between directional sensitivity and the overall event rate. In the 
event that directional detectors do indeed detect WIMP-induced recoil events, such studies 
will be crucial in forming a complete picture of the dark-matter particle and the structure 
of the Galactic halo. 

Acknowledgments 

S.K.L. thanks Marc Kamionkowski and Jens Chluba for useful comments. This research 
was supported at Caltech and Johns Hopkins University by DoE DE-FG03-92-ER40701. 
A. H. G. P. is supported by a Gary McCue Fellowship through the Center for Cosmology 
at UC Irvine and NASA Grant No. NNX09AD09G. 

A Event-Generation Procedure 

For simple velocity distributions, eq. (2.2) may be used to generate random scattering 
events and their corresponding nuclear-recoil momentum vectors directly from the analytic 
form of the recoil spectrum. However, for more general velocity distributions, an analytic 
form for the Radon transform may not be easily found. In such a case, it is then easier to use 
the velocity distribution to generate random incoming- WIMP velocity vectors, and then 
to randomly generate the corresponding nuclear-recoil momentum distribution by noting 



- 14 - 



that the elastic scattering is isotropic in the center-of-mass frame of the WIMP-nucleus 
system. 

We first note that the rate for scattering by a WIMP with velocity v is proportional 
to vf{v), as can be shown from eqs. (2.2) and (2.3). We thus draw randomly generated 
incoming- WIMP velocity vectors from the distribution proportional to vf{v). This is done 
using the method given in section 7.3 of ref. [84]. That is, we uniformly sample points in the 
4-dimensional space (v, (7), where the range of sampled velocity vectors v is determined by 
the truncation by the escape speed, and g S [0, max(t;/)] is at most the maximum value of 
f/(v). If a sampled point fails to satisfy the criteria vf{v) > g, then the sampled velocity 
vector V is discarded; otherwise, it is retained. The retained velocity vectors will then be 
distributed according to f/(v), as was desired. 

We then use these generated incoming- WIMP velocities to generate the corresponding 
nuclear recoil momenta. We first note that the recoil direction angular distribution is 
isotropic in the center-of-mass frame for non-relativistic elastic scattering events; however, 
this clearly implies that that the lab-frame recoil direction distribution is anisotropic, and 
will be a function of the incoming velocity v. To be specific, it is the distribution of the 
angle between v and q that is anisotropic in the lab frame. Thus, we could determine the 
appropriate distribution for each generated v, and then draw the recoil directions from 
these anisotropic distributions. However, it is easier to simply randomly draw the recoil 
direction from the isotropic center-of-mass frame distribution, and then transform it to the 
lab frame appropriately. 

Let the direction v of the incoming- WIMP lab-frame velocity vector be given by the 
spherical coordinates {6^, (j)^). Furthermore, let the angle between v and the nuclear recoil 
direction q be in the lab frame, and ^CM in the center-of-mass frame; these angles are 
related by 

cos0 = ^(1 + cos0cm)^/^. (A.l) 
v2 

Isotropic scattering in the center-of-mass frame implies that cos^cM is uniformly dis- 
tributed in the range [—1,1]; the distribution for cosO then follows. However, it also 
implies that the angle (j) between the plane defined by the lab-frame z-axis and v and the 
scattering plane (defined by q and v) is uniformly distributed in the range [0, 27r]. 

Thus, we determine q by first randomly generating the vector qo with spherical 
coordinates given by (^v + ^(^cm), </>¥)) which lies in the z — v plane and forms an angle 
with V. Rotating qo by the randomly generated angle (j) around the v axis then gives q 
with a distribution that is kinematically consistent with the incoming- WIMP velocity v. 
That is, q = R((/), v)qo, where in Cartesian coordinates 

(cos + ?;^(1 — cos 0) VxVy{l — cos (j)) — Vz sin</) ^^^^(l — cos (j)) + Vy sui(j)\ 
VxVy{l — cos (f) + Vz sin (f) cos (j) + Vy{l — cos (j)) VyVz{'^ — cos (f) — Vx sin cf) 
VxVz{^ — cos (j)) — Vy sin cf) VyVz{'^ — cos (f) + Vx sin cf) cos (j) + v1{l — cos (j)) j 

(A.2) 
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Finally, the nuclear-recoil momentum q = 2/icos^ is given by the kinematics. Thus, 
we have shown how to randomly generate the recoil momentum vectors q = g'q correspond- 
ing to randomly generated incoming- WIMP velocity vectors v drawn from a given velocity 
distribution /(v). 
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Table 2. Fiducial parameter values and flat prior ranges used to simulate data and perform likelihood analyses, respectively. Marginalized 
posterior probability modes and 68% minimum credible intervals are also given in a second set of rows for the halo-only 6-parameter, halo+stream, 
and halo+disk analyses. Parameters for the stream and disk components (denoted by a subscript C) are both given in the 5 leftmost columns. 




Energy (keV) 

Figure 1. Left: Simulated recoil spectrum for the halo-only 2-parameter and 3-parameter analyses 
in section 3.1.1, with binned signal events. The halo-only spectrum calculated from the fiducial 
values is also plotted, and is shown in solid black inside the energy sensitivity range. Right: 
Simulated recoil map for these analyses, in MoUweide projection. The black dot indicates the 
direction of the LSR, given in Galactic coordinates by (^lsr, ^lsr) = (90°, 0°). 
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Figure 2. Top row: Contour plots for the 2D posterior probability distribution in uiab-cn space, 
for the halo-only 2-parameter analyses with fixed in section 3.1.1. Analyses using energy- 
only, direction-only, and direction-t-energy information are shown. Red square markers indicate the 
fiducial values used in simulating the data. Black dots indicate the values used to generate the 
spectra and maps in figures 3 and 4. The marginalized posterior probability is shaded blue, with 
contours indicating 68% and 95% confidence levels. Bottom row: The same for the 3-parameter 
analyses assuming a flat mass prior. Note that the energy-only contours are larger in the (Th direction 
when only a flat prior on m-,^ is known instead of the exact value; in contrast, the direction-(-energy 
contours are relatively insensitive to lack of knowledge about m^. 



-27- 



^^lab = 50 km/s i^iab = 220 km/s fiab = 450 km/s 



0.500 



0.100 
0.050 

0.010 
0.005 



1.000 
0.500 



0.100 
0.050 

0.010 
0.005 



1.000 
0.500 



0.100 
0.050 

0.010 
0.005 



1.000 
0.500 



0.100 
0.050 



0.010 
0.005 



0.001 




■ 10 20 30 40 50 10 20 30 40 50 10 20 30 40 50 
Energy (keV) Energy (keV) Energy (keV) 

Figure 3. Recoil spectra showing the energy distribution (normalized to unity) of events in the 
5-50 keV energy range for a halo-only model, for the values of uiab and cth indicated by the black 
dots in figure 2. The similarity of some of the spectra here leads to the parameter degeneracy in 
the energy-only analysis shown in that figure. 
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Figure 4. Recoil maps showing the angular distribution (normalized to unity) of events in the 
5-50 keV energy range for a halo-only model, for the values of via,h and cth indicated in figure 2. 
Note the correspondence between the similarity of some of the maps and the parameter degeneracy 
in the direction-only analysis. Note also that as wiab increases, a deficit of events at the usual peak 
of the distribution appears, since the most energetic foward-scattcring events have energies greater 
than 50 keV at large viab- See also the discussion in rcf. [85]. 



- 29 - 




10 20 30 40 50 

Energy (keV) 

Figure 5. Simulated recoil spectrum for the halo-only 6-parameter analysis in section 3.1.2, with 
binned signal and background events. The fiducial halo-only spectrum with a flat background is 
also plotted. 




Figure 6. Simulated recoil maps for the halo-only 6-parameter analysis in section 3.1.2. Maps of 
the signal events due to WIMP-induced recoils in 5 different energy bins are shown, along with a 
map including both signal and background events. Note that the clustering of events tightens as 
their energy increases because of the kinematics of the scattering process. 
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Figure 7. Triangle plot showing ID and 2D posterior probability distributions over the full prior 
ranges, for the halo-only 6-paranieter analysis using only energy information in section 3.1.2. Red 
lines and square markers indicate the fiducial parameter values used in simulating the data. On the 
ID plots, 68% minimum credible intervals are shaded in blue. Parameter estimation is relatively 
poor. 
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Figure 8. The same as figure 7, but for the halo-only 6-parametcr analysis using direction-|-energy 
information in section 3.1.2. Parameter estimation is much improved over the energy-only analysis. 
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Figure 9. Simulated recoil spectrum for the halo+stream analysis in section 3.2. 
lialo+stream spectrum with a flat background is also plotted. 
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Figure 10. Simulated recoil maps for the halo+stream analysis in section 3.2, as in figure 6. 
Black dots indicate the direction of the LSR and the stream. The cluster of stream events is most 
easily seen in the top-right panel, although it is slightly visible as a wider ring-like feature at lower 
energies. 
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Figure 11. Triangle plot for the halo+stream analysis in section 3.2, using dircction+energy information. Parameter estimation 
is relatively good. 
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Figure 12. Simulated recoil spectrum for the halo+disk analyses in section 3.3. The fiducial 
halo+disk spectrum with a flat background is also plotted. The shaded light-blue bars indicate the 
binned events below the 5-keV energy threshold included in the second 0-50 keV analysis. 






Figure 13. Simulated recoil maps for the halo+disk analyses in section 3.3, as in figure 6. Black 
dots indicate the common direction of the LSR and the disk. Events over the full 0-50 keV energy 
range are shown. 
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Figure 14. Triangle plot for the 0-50 keV halo+disk analysis in section 3.3. Parameter estimation is fair, but some of the disk 
parameters are recovered with less accuracy. 
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Figure 15. Triangle plot for the 5-50 keV halo+disk analysis in section 3.3. The disk parameters are poorly estimated, 
demonstrating the need for improved energy thresholds. 



